home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / CUGUK / APPLICAT / C050.ZIP / SFSSRC.ZIP / SFS.SRC / AS / AS_SPJ.C < prev    next >
C/C++ Source or Header  |  1991-11-04  |  14KB  |  588 lines
  1. /***************************************************************
  2.  
  3.     as_spj.c        Spherical projection routines
  4.             for astronomy (as) subsystem
  5.  
  6.             Copyright (c) 1991, Ted A. Campbell
  7.  
  8.             Bywater Software
  9.             P. O. Box 4023
  10.             Duke Station
  11.             Durham, NC  27706
  12.  
  13.             email: tcamp@hercules.acpub.duke.edu
  14.  
  15.     Copyright and Permissions Information:
  16.  
  17.     All U.S. and international copyrights are claimed by the
  18.     author. The author grants permission to use this code
  19.     and software based on it under the following conditions:
  20.     (a) in general, the code and software based upon it may be
  21.     used by individuals and by non-profit organizations; (b) it
  22.     may also be utilized by governmental agencies in any country,
  23.     with the exception of military agencies; (c) the code and/or
  24.     software based upon it may not be sold for a profit without
  25.     an explicit and specific permission from the author, except
  26.     that a minimal fee may be charged for media on which it is
  27.     copied, and for copying and handling; (d) the code must be
  28.     distributed in the form in which it has been released by the
  29.     author; and (e) the code and software based upon it may not
  30.     be used for illegal activities.
  31.  
  32. ***************************************************************/
  33.  
  34. #include "stdio.h"
  35. #include "math.h"
  36. #include "bw.h"
  37. #include "as.h"
  38.  
  39. #ifdef __STDC__
  40. #include "malloc.h"
  41. #else
  42. extern char * malloc();
  43. #endif
  44.  
  45. /***************************************************************
  46.  
  47.     spj_readspd()            Read an SPJ datafile
  48.  
  49. ***************************************************************/
  50.  
  51. spj_readspd( datafile, start, end )
  52.    char *datafile;
  53.    struct spj_pt *start, *end;
  54.    {
  55.    FILE *data;
  56.    register int c;
  57.    static char read_buffer[ 128 ];
  58.    static int co;
  59.    static double la, lo, di;
  60.    struct spj_pt *current, *new;
  61.  
  62.    /* Open the datafile */
  63.  
  64.    if ( ( data = fopen( datafile, "rb" )) == NULL )
  65.       {
  66.       sprintf( bw_ebuf, "Failed to open datafile %s", datafile );
  67.       bw_error( bw_ebuf );
  68.       return BW_ERROR;
  69.       }
  70.    else
  71.       {
  72.       sprintf( bw_ebuf, "Reading spherical projection datafile %s",
  73.      datafile );
  74.       bw_message( bw_ebuf );
  75.       }
  76.  
  77.    /* Get storage for the first point */
  78.  
  79.    if ( ( current = ( struct spj_pt  *) malloc( sizeof( struct spj_pt ) )) == NULL )
  80.       {
  81.       sprintf( bw_ebuf, "Failed to allocate memory reading %s", datafile );
  82.       bw_error( bw_ebuf );
  83.       return BW_ERROR;
  84.       }
  85. #ifdef  OLD_DEBUG
  86.    else
  87.       {
  88.       sprintf( bw_ebuf, "Memory for %s allocated ok",
  89.      datafile );
  90.       bw_message( bw_ebuf );
  91.       }
  92. #endif
  93.  
  94.     start->next = current;
  95.  
  96.    /* Read the data from the file */
  97.  
  98.    c = 0;
  99.    while( feof( data ) == 0 )
  100.       {
  101.  
  102.       as_fgets( read_buffer, 127, data );
  103.       sscanf( read_buffer, "%d%lf%lf%lf", &co, &la, &lo, &di );
  104.       current->code = co;
  105.       current->latitude  = la;
  106.       current->longitude = lo;
  107.       current->radius    = di;
  108.  
  109. #ifdef  OLD_DEBUG
  110.       sprintf( bw_ebuf, "[pr:] spj_readsd(): %d: %.1lf %.1lf",
  111.       c, la, lo );
  112.       bw_message( bw_ebuf );
  113. #endif
  114.  
  115.       /* Now get a new 'current' */
  116.  
  117.       if ( ( new = (struct spj_pt  *)
  118.      malloc( sizeof( struct spj_pt ) )) == 0 )
  119.      {
  120.      sprintf( bw_ebuf, "Failed to allocate memory reading %s", datafile );
  121.      bw_error( bw_ebuf );
  122.      return BW_ERROR;
  123.      }
  124. #ifdef  OLD_DEBUG
  125.       else
  126.      {
  127.      sprintf( bw_ebuf, "Memory for %s allocated ok",
  128.         datafile );
  129.      bw_message( bw_ebuf );
  130.      }
  131. #endif
  132.  
  133.       /* Reassign linkage from current to new */
  134.  
  135.       if ( feof( data ) == FALSE )
  136.      {
  137.      current->next = new;
  138.      current       = new;
  139.      current->next = end;
  140.      }
  141.  
  142.       if ( ferror( data ) != FALSE )
  143.      {
  144.      sprintf( bw_ebuf, "Failed reading %s", datafile );
  145.      bw_error( bw_ebuf );
  146.      return BW_ERROR;
  147.      }
  148.       ++c;
  149.       }
  150.  
  151.    fclose( data );
  152.    return TRUE;
  153.    }
  154.  
  155. /***************************************************************
  156.  
  157.     spj_calc()      Calculate a linked set of SPJ
  158.             point structures
  159.  
  160. ***************************************************************/
  161.  
  162. spj_calc( start, end, vlat, vlon, vdis, forced_radius, rotation, mode, poll )
  163.    struct spj_pt *start, *end;
  164.    double vlat, vlon, vdis;
  165.    double forced_radius;
  166.    double rotation;
  167.    int mode;
  168.    int (*poll)();
  169.    {
  170.    struct spj_pt *current;
  171.    static double x, y;
  172.    static int side;
  173.    register int d;
  174.    double radius;
  175.  
  176. #ifdef  DEBUG
  177.    if ( ( vlat < -90.0 ) || ( vlat > 90.0 ))
  178.       {
  179.       sprintf( bw_ebuf, "[pr:] spj_calc() received vlat = %lf", vlat );
  180.       bw_error( bw_ebuf );
  181.       return;
  182.       }
  183.    if ( ( vlon < -180.0 ) || ( vlon > 180.0 ))
  184.       {
  185.       sprintf( bw_ebuf, "[pr:] spj_calc() received vlon = %lf", vlon );
  186.       bw_error( bw_ebuf );
  187.       return;
  188.       }
  189.    if ( vdis < 0.0 )
  190.       {
  191.       sprintf( bw_ebuf, "[pr:] spj_calc() received vdis = %lf", vdis );
  192.       bw_error( bw_ebuf );
  193.       return;
  194.       }
  195.    if ( ( rotation < 0.0 ) || ( rotation > 360.0 ))
  196.       {
  197.       sprintf( bw_ebuf, "[pr:] spj_calc() received rotation = %lf", rotation );
  198.       bw_error( bw_ebuf );
  199.       return;
  200.       }
  201.    if ( ( mode < 0 ) || ( mode > 2 ))
  202.       {
  203.       sprintf( bw_ebuf, "[pr: spj_calc() received mode = %d", mode );
  204.       bw_error( bw_ebuf );
  205.       return;
  206.       }
  207. #endif
  208.  
  209.    d = 0;
  210.    current = start->next;
  211.    while( current != end )
  212.       {
  213.  
  214.       ++d;
  215.  
  216.       if ( forced_radius == 0.0 )
  217.      {
  218.      radius = current->radius;
  219.      }
  220.       else
  221.      {
  222.      radius = forced_radius;
  223.      }
  224.  
  225. #ifdef    OLD_DEBUG
  226.       sprintf( bw_ebuf, " calculate %d: lat %0.2lf, lon %0.2lf, rad %0.2lf      >>",
  227.      d, current->latitude, current->longitude, radius );
  228.       bw_message( bw_ebuf );
  229.       kb_rx();
  230. #endif
  231.  
  232.       spj_point( current->latitude, current->longitude,
  233.      radius, vlat, vlon, vdis, rotation,
  234.      mode, &x, &y, &side );
  235.       current->x    = x;
  236.       current->y    = y;
  237.       current->side = side;
  238.  
  239. #ifdef    OLD_DEBUG
  240.       sprintf( bw_ebuf, " x = %0.2lf, y = %0.2lf, v = %d      >>",
  241.      current->x * ACC_FAC,
  242.      current->y * ACC_FAC,
  243.      current->side );
  244.       bw_message( bw_ebuf );
  245.       kb_rx();
  246. #endif
  247.  
  248.       /* call poll function */
  249.  
  250.       (*poll) ();
  251.  
  252.       current = current->next;
  253.       }
  254.    }
  255.  
  256. /***************************************************************
  257.  
  258.     spj_point()     Calculate a spherical projection point
  259.  
  260.     This routine is, in a sense, the heart of the Space
  261.     Flight Simulator and other spherical-projection programs
  262.     related to it, such as the Space Flight Atlas.
  263.  
  264. ***************************************************************/
  265.  
  266. spj_point( plat, plon, prad, vlat, vlon, vdis, rotation, mode, x, y, side )
  267.    double plat;       /* latitude of point in object coordinates */
  268.    double plon;       /* longitude of point in object coordinates */
  269.    double prad;       /* "radius" i.e., distance from center of
  270.            object, of the point */
  271.    double vlat;       /* latitude of viewer */
  272.    double vlon;       /* longitude of viewer */
  273.    double vdis;       /* distance of viewer from center of object */
  274.    double rotation;   /* rotation of object in degrees, counter-
  275.            clockwise;
  276.    int mode;          /* mode of calculation, i.e., side(s) to calculate */
  277.    double *x;         /* return x viewer coordinate in degrees of arc */
  278.    double *y;         /* return y viewer coordinate in degrees of arc */
  279.    int *side;         /* return side of sphere */
  280.    {
  281.    static double old_vlat = 361.0;
  282.    static double old_vlon = 361.0;
  283.    static double old_vdis = -1.0;
  284.    static double old_plat = 361.0;
  285.    static double old_plon = 361.0;
  286.    static double old_prad = -1.0;
  287.    static double old_rot  = 361.0;
  288.    static double x_lon    = 361.0;
  289.    static double x_rsin, x_rcos;
  290.    static double x_cosplat, x_sinplat;
  291.    static double x_cosvlat, x_sinvlat;
  292.    static double x_angrad;
  293.    double x_cosc, x_cosl;
  294.    double x_temp, y_temp;
  295.  
  296. #ifdef  OLD_DEBUG
  297.    sprintf( bw_ebuf, "[pr:] spj_point() received: plat %0.2lf, plon %0.2lf, prad %0.2lf      >>",
  298.       plat, plon, prad );
  299.    bw_message( bw_ebuf );
  300.    kb_rx();
  301.    sprintf( bw_ebuf, "[pr:] spj_point() received: vlat %0.2lf, vlon %0.2lf, vdis %0.2lf      >>",
  302.       vlat, vlon, vdis );
  303.    bw_message( bw_ebuf );
  304.    kb_rx();
  305.    sprintf( bw_ebuf, "[pr:] spj_point() received: rot %0.2lf, mode %d      >>",
  306.       rotation, mode );
  307.    bw_message( bw_ebuf );
  308.    kb_rx();
  309. #endif
  310.  
  311. #ifdef  DEBUG
  312.    if ( ( vlat < -90.0 ) || ( vlat > 90.0 ))
  313.       {
  314.       sprintf( bw_ebuf, "[pr:] spj_calc() received vlat = %lf", vlat );
  315.       bw_error( bw_ebuf );
  316.       return;
  317.       }
  318.    if ( ( vlon < -180.0 ) || ( vlon > 180.0 ))
  319.       {
  320.       sprintf( bw_ebuf, "[pr:] spj_calc() received vlon = %lf", vlon );
  321.       bw_error( bw_ebuf );
  322.       return;
  323.       }
  324.    if ( vdis < 0.0 )
  325.       {
  326.       sprintf( bw_ebuf, "[pr:] spj_calc() received vdis = %lf", vdis );
  327.       bw_error( bw_ebuf );
  328.       return;
  329.       }
  330.    if ( ( plat < -90.0 ) || ( plat > 90.0 ))
  331.       {
  332.       sprintf( bw_ebuf, "[pr:] spj_calc() received plat = %lf", plat );
  333.       bw_error( bw_ebuf );
  334.       return;
  335.       }
  336.    if ( ( plon < -180.0 ) || ( plon > 180.0 ))
  337.       {
  338.       sprintf( bw_ebuf, "[pr:] spj_calc() received plon = %lf", plon );
  339.       bw_error( bw_ebuf );
  340.       return;
  341.       }
  342.    if ( prad < 0.0 )
  343.       {
  344.       sprintf( bw_ebuf, "[pr:] spj_calc() received prad = %lf", prad );
  345.       bw_error( bw_ebuf );
  346.       return;
  347.       }
  348.    if ( ( rotation < 0.0 ) || ( rotation > 360.0 ))
  349.       {
  350.       sprintf( bw_ebuf, "[pr:] spj_calc() received rotation = %lf", rotation );
  351.       bw_error( bw_ebuf );
  352.       return;
  353.       }
  354.    if ( ( mode < 0 ) || ( mode > 2 ))
  355.       {
  356.       sprintf( bw_ebuf, "[pr: spj_calc() received mode = %d", mode );
  357.       bw_error( bw_ebuf );
  358.       return;
  359.       }
  360. #endif
  361.  
  362.    /* if longitude is changed, calculate its sine and cosine */
  363.  
  364.    if ( plon != old_plon )
  365.       {
  366.       x_lon    = spj_meridian( vlon, plon );
  367.       old_plon = plon;
  368. #ifdef  OLD_DEBUG
  369.       bw_message( "[pr:] spj_point() recalculated plon" );
  370.       kb_rx();
  371. #endif
  372.       }
  373.  
  374.    /* if viewer latitude is changed, calculate its sine and cosine */
  375.  
  376.    if ( vlat != old_vlat )
  377.       {
  378.       x_cosvlat = vpt_cos( vlat );
  379.       x_sinvlat = vpt_sin( vlat );
  380.       old_vlat = vlat;
  381. #ifdef  OLD_DEBUG
  382.       bw_debug( "[pr:] spj_point() recalculated vlat" );
  383. #endif
  384.       }
  385.  
  386.    /* if point latitude is changed, calculate its sine and cosine */
  387.  
  388.    if ( plat != old_plat )
  389.       {
  390.       x_cosplat = vpt_cos( plat );
  391.       x_sinplat = vpt_sin( plat );
  392.       old_plat = plat;
  393. #ifdef  OLD_DEBUG
  394.       bw_debug( "[pr:] spj_point() recalculated plat" );
  395. #endif
  396.       }
  397.  
  398.    /* calculations to determine side */
  399.  
  400.    x_cosl    = vpt_cos( x_lon ) * x_cosplat;
  401.    x_cosc    = x_sinvlat * x_sinplat + x_cosvlat * x_cosl;
  402.  
  403.    /* Note side of sphere, and return if this side is not
  404.       to be calculated now */
  405.  
  406.    if ( x_cosc < 0 )
  407.       {
  408.       *side = SPJ_FARSIDE;
  409.       if ( mode == SPJ_NEARSIDE )
  410.      {
  411.      return;
  412.      }
  413.       }
  414.    else
  415.       {
  416.       *side = SPJ_NEARSIDE;
  417.       if ( mode == SPJ_FARSIDE )
  418.      {
  419.      return;
  420.      }
  421.       }
  422.  
  423.    /* if point radius or viewer distance is changed, recalculate
  424.       the angular radius */
  425.  
  426.    if ( ( prad != old_prad ) || ( vdis != old_vdis ))
  427.       {
  428.       x_angrad = spj_angrad( vdis, prad );
  429.       old_prad = prad;
  430.       old_vdis = vdis;
  431. #ifdef  OLD_DEBUG
  432.       bw_message( "[pr:] spj_point() recalculated angular radius" );
  433.       kb_rx();
  434. #endif
  435.       }
  436.  
  437.    /* Assign unrotated values */
  438.  
  439.    *x = x_angrad * ( x_cosplat * vpt_sin( x_lon ) );
  440.    *y = x_angrad * ( x_cosvlat * x_sinplat - x_sinvlat * x_cosl );
  441.  
  442.    /* if rotation has changed, recalculate its sine and cosine */
  443.  
  444.    if ( rotation != old_rot )
  445.       {
  446.       x_rsin  = vpt_sin( rotation );
  447.       x_rcos  = vpt_cos( rotation );
  448.       old_rot = rotation;
  449. #ifdef  OLD_DEBUG
  450.       bw_message( "[pr:] spj_point() recalculated rotation" );
  451.       kb_rx();
  452. #endif
  453.       }
  454.  
  455.    /* Rotate the image */
  456.  
  457.    x_temp = *x;
  458.    y_temp = *y;
  459.    *x = 0.0 - (( x_temp * x_rcos ) - ( y_temp * x_rsin ));
  460.    *y = ( x_temp * x_rsin ) + ( y_temp * x_rcos );
  461.  
  462. #ifdef  OLD_DEBUG
  463.    sprintf( bw_ebuf, "[pr:] spj_point() returns: x %0.2lf, y %0.2lf, side %d      >>",
  464.       *x, *y, *side );
  465.    bw_message( bw_ebuf );
  466.    kb_rx();
  467. #endif
  468.  
  469.    }
  470.  
  471. /***************************************************************
  472.  
  473.     spj_meridian()  Calculate the distance in degrees
  474.             between a point's longitude and the
  475.             viewer's longitude
  476.  
  477. ***************************************************************/
  478.  
  479. double
  480. spj_meridian( vlon, plon )
  481.    double vlon;       /* viewer's longitude */
  482.    double plon;       /* point longitude */
  483.    {
  484.    double retval;
  485.  
  486.    retval = vlon - plon;
  487.    if ( retval < -180 )
  488.       {
  489.       retval += 360;
  490.       }
  491.    if ( retval > 180 )
  492.       {
  493.       retval -= 360;
  494.       }
  495.    return retval;
  496.    }
  497.  
  498. /***************************************************************
  499.  
  500.    spj_angrad()    Calculate the angular radius of an object
  501.      given the radius of the object and the
  502.      distance to the viewer
  503.  
  504. ***************************************************************/
  505.  
  506. double
  507. spj_angrad( distance, radius )
  508.    double distance, radius;
  509.    {
  510.  
  511. #ifdef  DEBUG
  512.    if ( distance == 0.0 )
  513.       {
  514.       sprintf( bw_ebuf, "[pr:] spj_angrad() received distance = %lf",
  515.      distance );
  516.       bw_error( bw_ebuf );
  517.       return 1.0;
  518.       }
  519. #endif
  520.  
  521.    return ( RAD_DEG * atan( radius / distance ));
  522.    }
  523.  
  524. /***************************************************************
  525.  
  526.    spj_degfix()    Fix degrees, i.e., 0 < d < 360
  527.  
  528. ***************************************************************/
  529.  
  530. double
  531. spj_degfix( n )
  532.    double n;
  533.    {
  534.    double retval;
  535.  
  536.    retval = n;
  537.    while ( retval < 0 )
  538.       {
  539.       retval += 360;
  540.       }
  541.    while ( retval > 360 )
  542.       {
  543.       retval -= 360;
  544.       }
  545.  
  546. #ifdef  OLD_DEBUG
  547.    sprintf( bw_ebuf, "spj_degfix(): rec %.2lf, ret %.2lf", n, retval );
  548.    bw_debug( bw_ebuf );
  549. #endif
  550.  
  551.    return retval;
  552.    }
  553.  
  554. /***************************************************************
  555.  
  556.     as_fgets()      Gets a string from a bitstream,
  557.             excluding strings that begin with ';'
  558.  
  559.     This function acts exactly as the standard Unix/C
  560.     fgets() function, except that it rejects lines that
  561.     begin with a semicolon.  This alllows us to imbed
  562.     comment lines in as datafiles by beginning them with
  563.     a semicolon.
  564.  
  565. ***************************************************************/
  566.  
  567. char *
  568. as_fgets( s, n, stream )
  569.    char *s;
  570.    int n;
  571.    FILE *stream;
  572.    {
  573.    char *r;
  574.  
  575.    r = (char *) 1;
  576.    s[ 0 ] = ';';
  577.    while( s[ 0 ] == ';' )
  578.       {
  579.       r = fgets( s, n, stream );
  580.       if ( r == NULL )
  581.      {
  582.      return NULL;
  583.      }
  584.       }
  585.    return r;
  586.    }
  587.  
  588.